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^ ■ Abstract 
o 

c/3 ■ We apply first principles computational techniques to analyze the two-electron, multi-step, elec- 

[ trochemical reduction of CO2 to CO in water using cobalt porphyrin as a catalyst. Density Func- 



tional Theory calculations with hybrid functionals and dielectric continuum solvation are used to 



^ ' determine the steps at which electrons are added. This information is corroborated with ah initio 

^ : 

10 . molecular dynamics simulations in an explicit aqueous environment which reveal the critical role 

o 

of water in stabilizing a key intermediate formed by CO2 bound to cobalt. Using potential of mean 
force calculations, the intermediate is found to spontaneously accept a proton to form a carboxylate 
acid group at pH<9.0, and the subsequent cleavage of a C-OH bond to form CO is exothermic and 
associated with a small free energy barrier. These predictions suggest that the proposed reaction 

XI 

%-{ ' mechanism is viable if electron transfer to the catalyst is sufficiently fast. The variation in cobalt 

ion charge and spin states during bond breaking, DFT-I-U treatment of cobalt 3d orbitals, and the 
need for computing electrochemical potentials are emphasized. 



I. INTRODUCTION 



CO2 capture from flue gas and its conversion to useful products, including fuel molecules, 
has emerged as an important paradigm for a carbon- neutral economy. At discussed in the 
preceeding paper of this series (henceforth "paper I"),— high (~70%) yield of carbon monox- 
ide (CO) has been demonstrated in cobalt macro cycle- catalyzed electrochemical reduction 
of carbon dioxide (CO2) in water at applied voltage of about -1.0 volt.-"— The mechanisms 
of CO2 reduction in non-aqueous solvents, for which much more negative potentials are 
needed, have been examined using a variety of methods.—"— 

Co(I)P-catalyzed CO2 reduction in water and protic solvents, which requires a much 
less negative voltage for the onset of reaction than in organic solvent,-"— has received less 
fundamental studies. The present theoretical work focuses on the mechanism of this elec- 
trochemical reaction in aqueous media. As discussed in the preceeding paper in this series22 
(henceforth "Paper I" ) which examines the structures, energetics, and charge states of reac- 
tion intermediates in detail, the reaction likely takes place in the following logical sequence 
of steps: 

C0P + CO2 ^ C0PCO2; (1) 
C0PCO2 + H+ ^ CoPCOOH; (2) 
CoPCOOH ^ CoPCO + OH- ; (3) 
CoPCO ^ CoP + CO. (4) 

"CoP" will henceforth denote cobalt porphine, the simplest porphyrin species, adopted in 
this work for ease of calculations.— Key questions to be addressed include: (1) why the 
reaction, which involves protonation to form carboxylate acid motifs (COOH) with pi^a 
typically on the order of 4.5, readily proceeds despite the fact that pH> 7 in experiments; 
(2) whether all steps are thermodynamically downhill; (3) whether the free energy barriers 
of the intermediate steps are low enough to be consistent with the observed reaction rate; 
and (4) at what stages the two electrons are added. 

In Eqs. [T]|l] above, we have intentionally left out the charge states of the intermediates 
as yet unassigned in experiments. Since the two electrons can be added at any step(s), 
numerous mechanisms are consistent with these equations. The sequence of electron injection 
is governed by the redox potentials ($redox) of the pertinent reaction intermediates relative 



to the applied voltage. The half-cell potentials of various charge states of cobalt porphyrins 
are known/^ but not those for the CO- and CO2- bound complexes. 

In this work, we use a combination of ab initio molecular dynamics (AIMD)^ techniques 
and Density Functional Theory (DFT) calculations with various exchange-correlation func- 
tional and the polarizable continuum model^^ (hereafter referred to as "DFT+pcm") to 
study Eqs. [T]|H These methods inform and support each other. The more economical 
DFT+pcm calculations approximate the aqueous environment as a dielectric continuum,— 
allowing a global overview of the entire electrochemical reaction and a survey of the numer- 
ous reaction intermediates. Thus, DFT+pcm results reported in Paper I are used to ex- 
tract $rcdox for all possible intermediates. Redox potentials consist of hydration free energy 
(AGhyd) and ionization potential/electron affinity contributions. While $redox and AGhyd 
have also been calculated with the AIMD method and explicit treatment of the aqueous 
environment,—"— these have so far been limited to monoatomic ions and molecules much 
smaller than porphyrins comprising a minimum of 37 atoms.^'' Our B3LYP DFT+pcm re- 
dox potentials predict that CoPCOOH" is the key intermediate. To corroborate this $redox 
conclusion, we also conduct more costly AIMD simulations with explicit water molecules. 
The hydration structures of the reaction intermediates predicted therein help explain the 
redox potential trends. 

Next, we perform AIMD calculations on the individual reactions identified as the key 
steps in DFT+pcm calculations. Although it is possible that the electron addition steps of 
the two-electron reduction of CO2 are rate-limiting,—*^ modeling the electron transfer rate 
via Marcus theory approach^^ may not answer the most interesting scientific questions. This 
is because the electron injection rate almost certainly depends on engineering aspects such 
as the electrical contact between the gas diffusion electrode and the polymerized catalyst. 
Besides, it is necessary to demonstrate the viability of other steps in which electron transfer 
is not involved. We omit the electrode, focus on cobalt porphine molecule dispersed in 
water, and use AIMD to study two reactions involving [Co(II)PCOOH]~: deprotonation 
(the reverse of Eq. |2]) and the cleavage of the C-OH bond (Eq. |3]). The protonation reaction 
proves critical to the efficient removal of one of the CO2 oxygen atoms and reduction of the 
carbon atom from the +4 to the +2 formal charge state. The pK^, of [Co(II)PCOOII]^ will 
be determined using AIMD potential of mean force (PMF) in a manner previously applied 
to silanol groups on silica surfaces. As for the free energy change and barrier associated 



with the breaking of the C-OH bond (Eq. [3]) prior to releasing CO (Eq. |4]), there is as yet no 
experimental reaction rate for direct comparison. The turnover rate per catalytic active site 
is obscured by the undetermined proportion of active CoP molecules actually participating 
in the reactions. But our PMF with an approximate reaction coordinate yields a barrier 
which suggests that CO gas evolution should proceed readily if electron transfer from the 
electrode is fast. It is known experimentally that Co(II)P binds weakly to CO;^"^ our gas 
phase calculations also shows that the Co(l)P-CO complex is weakly bound. Therefore once 
the C-OH bond is severed to form a 0H~, the subsequent steps (Eq. H]) should be exothermic, 
fast, and non-rate-determining, and they do not require further theoretical studies. 

This work provides fundamental understanding specific to the CO2 reduction mechanism 
in water using macromolecule catalysts, highlights the critical role played by the protic 
solvent water in lowering the voltage needed for the reaction, and may shed light on ways 
to further improve and modify the cobalt porphyrin catalyst. It is also of general interest to 
the fledging field of computational electrochemistry. The multistep nature of the reaction 
emphasizes the importance of accurate calculation of redox potentials.—""— The protonation 
reaction (Eq. [2]) is accompanied with a change in the Co ion charge state, which can be 
considered a form of coupled proton/electron process.— 1^ We show that this can induce 
hysteresis in AIMD simulations. Finally, this study in explicit water, comprising AIMD 
trajectories exceeding 200 ps in duration, is greatly facilitated by the use of an empirical 
DFT+U method^ that can treat the localized 3(i-orbitals of the Co ion accurately without 
resorting to more costly theoretical methods such as hybrid functionals. 

This paper is organized as follows: The theoretical methods used are discussed in Sec. 2. 
Section 3 describes the redox potential predictions which determine the electron addition 
steps, and then focuses on the deprotonation and C-0 bond cleavage reactions of the key 
intermediate [Co(II)PCOOH]^. Section 4 compares the method used in the present work to 
our previous C02-related simulations;— it further looks to the future and briefly discusses 
new computational techniques that may facilitate the modeling of demanding electrochem- 
ical processes. The accuracy of electronic structure calculations which underlie our mecha- 
nistic predictions will also be addressed.— Section 5 concludes and summarizes the paper. 
An appendix discusses minor hysteresis issues encountered in one part of the calculations. 



II. METHOD 



All B3LYP or PBE plus dielectric continuum calculations apply the Gaussian suite of 
programs version g03.^° All DFT+U calculations use the VASP code.— The Supporting 
Information document (SI) Sec. SI provides a brief comparison between these packages. 

DFT plus dielectric continuum (DFT+pcm) calculations apply the Becke-3-parameter- 
Lee- Yang-Parr (B3LYP)^i^ or the Perdew-Burke-Ernzerhof (PBE)^^ functional, and the 
PCM dielectric continuum model.— Other details are described in Paper I.^^ Zero-point 
energy (ZPE) contributions to Eqs. [2] and [3] are estimated at T=0 K, also using Gaussian^ 
and the 6-31-l-G* basis set. The H20-0H~ complex is assumed to be the H"*" accepting or 
0H~ containing species in ZPE calculations. 

The redox potential of the A'^^+^^^/A"^ couple is the electron affinity (EA) of A"~ plus the 
difference in hydration free energies (AGhyd difference between A("+i)- and A""). Both EA 
and AGhyd are readily obtained from DFT-|-pcm total free energies that include zero point 
energy corrections and finite temperature contributions. Where available, DFT+U $rcdox 
are estimated by substituting ground state DFT+U energies for B3LYP ones but retaining 
B3LYP dielectric hydration and ZPE information.— All reported $redox are referenced to 
the accepted value of 4.44 volt for the standard hydrogen electrode half-cell potential. The 
effect of the choice of DFT functionals on $redox will be addressed in Sec. [TVl 

Spin-polarized AIMD simulations apply Vienna ab-initio simulation program (VASP),— 
the PBE functional, F-point Brillouin zone sampling, 400 eV planewave energy cutoff, deu- 
terium masses for all protons to allow Born-Oppenheimer dynamics time steps of 0.25 fs, 
a 10~^ eV energy convergence criterion, and T=425 K NVT conditions using a Nose ther- 
mostat. At T=400 K, the PBE functional yields average pure water structure consistent 
with experimentally observed water g{r) at T=300 K.^^ Since the porphine ring exhibits 
significant ruffling fluctuations that explore a large conflgurational space, we raise the tem- 
perature by an extra 25 K to promote better sampling statistics. AIMD simulations apply 
13.64 Ax 13.64 Ax 13.64 A simulation cells which contain a CoPCOOH complex and 71 H2O 
molecules.— 

Semi-local functionals such as PBE are generally inadequate for treating flrst row transi- 
tion metal complexes like cobalt,—"^ although they are more successful with transition metal 
dimers.— In general, transition metal complexes have been a challenge to DFT methods.^- 



We apply the following reaction to benchmark the prediction energetics: 



[Co(III)P]+ + CO ^ [Co(III)P - CO] 



(5) 



While Co(II)P itself seems well represented by DFT functionals,— Paper I has shown that 
the widely used PBE and B3LYP functionals predict Eq. [5]binding energies which differ from 
the experimental valued by about 0.5 eV in opposite directions. To deal with this problem, 
we augment the PBE functional with DFT+U— applied to the partially occupied 3(i-orbitals 
of the Co ion. With a judicial choice of the U parameter, this approach has been shown 
to give accurate predictions for organometallic compounds, although agreement between 
theory and experiments has not been universal.—"— Using VASP PAW pseudopotentials,— 
setting f/=2.5 eV yields a binding energy for Eq. [5] that agrees with experimental data 
reported for the Co(III)TPP-CO complex.— i^i^ At T=0 K, these two complexes exhibit 
very similar VASP/PBE binding energies of 1.319 and 1.351 eV respectively; they differ by 
only 0.74 kcal/mol, suggesting that the experimental Co(III)TPP-CO binding free energy is 
a good metric for benchmarking the theoretical Co(III)P-CO predictions. As neither CO2 
nor CO strongly binds to cobalt porphyrins at most accessible Co charge states, Eq. |5]is the 
only binding constant available in the experimental literature as a benchmark. While this 
value of U may not be optimal for all Co charge states, referencing to Eq. |5] appears the 
most justifiable empirical route. The accuracy of AIMD simulations depend on the DFT+U 
method used, as will be discussed in Sec. HVl 

As reported in Paper I,— the optimal spin states of all gas phase B3LYP and PBE CoP 
complexes are low spin except [Co(I)PCOOH]^~ and the unligated [Co(III)P] + , both of 
which are triplets. DFT+U calculations at T=0 K yield the same optimal spin states. 

Equation [2] involves calculating the pi^a of [Co(II)PCOOH]". pKg, = — log^g 
exp(— Z^AG'^'')) has been successfully computed for molecules and surfaces in liquid water 
using the AIMD technique.— 1^"— Here (3 is l/ksT, and AG^^^ is the standard state depro- 
tonation free energy. 



Co denotes 1.0 M concentration, R is the reaction coordinate, A{R) is a phase space factor to 
be discussed in Sec. lIIIl i?cut is the cutoff distance delimiting the reaction and product valleys 
in the free energy landscape, and W{R) is the potential of mean force (PMF), referenced 





such that W{R) = as i? — )• oo. -Rcut is taken as the onset of the plateau where W{R) — )■ 0. 
The umbrella sampling technique^ and a 4-atom reaction coordinate 

R = Ri — R2 — R3 (7) 

are applied to compute W{R). Here Ri is the distance between the COOH acid proton and 
the oxygen atom on the designated proton-accepting H2O, while R2 and -R3 are the dis- 
tances between this O atom and the two protons originally on the designated H2O molecule, 
respectively.— As R2 and -R3 are about 1.0 A for intact 0-H covalent bonds, it can be read- 
ily inferred that i? ~ — 1.0 A is consistent with deprotonation and CO^/HsO^ contact ion 
pair formation while R > —0.4 A indicates an intact CO-H bond.'^^ Designating a special 
H2O molecule can be done without loss of generality because all water molecules are inter- 
changeable and only one is at any time is close enough to the acid proton to be considered 
a potential proton acceptor. Harmonic potentials of the form 

U{R) = B{R-Ro)\ (8) 

with B values between 2 to 4 eV/A^ are applied to 7 umbrella windows with Ro spanning 
the range of R to be sampled. We reference the pJ^a of [CoPCOOH]^ relative to the 
free energy of water-autoionization,''^ computed using the same reaction coordinate and 
elevated temperature and assumed to exhibit pK^=lA. As an validation test, our AIMD 
pKs, methodology has been applied to formic acid in water, yielding an acidity constant 
within a fraction of a pH unit of the experimental value (SI Sec. S2). 

A "reflecting boundary condition" potential V{Ro}i) sets an approximate 1.2 A distance 
of closest approach between water protons and the hydroxyl oxygen atom. It preserves 
the identity of the deprotonated -COOH (or H2O in the case of water auto-ionization) by 
preventing proton transfer via the Grotthuss mechanism. V^(-Roh) = B{Ron — RiY, where 
5=200 eV/A"^ and -Ri=1.3 A, is imposed whenever i?oH < Ri- Ron is the distance between 
the hydroxyl (COH) oxygen and all H"*" other than the original COOH proton. Related 
boundary potentials have been applied to AIMD simulations of other chemical reactons.— 
V{Ron) is only necessary in the deprotonation umbrella sampling window with the most 
negative R. As this leftmost umbrella sampling window exhibits a W{R) variation of only 
~ 0.6 kcal/mol (see Sec. IIIip . the effect of V{Ro}i) should be small — much less than 
0.6 kcal/mol. 



FIG. 1: Snapshot of (a) [Co(I)P] ; (b) [Co(I)P]^ (see footnote l9ll for description of charge states; 
(c) [Co(I)PC02]"; (d) [Co(I)PC02]2-; (e) Co(II)PCOOH; (f) [Co(II)PCOOH]-; in water. Panels 
(a)-(b) show the absence of hydrogen bonding between the negatively charged CoP and water. 
In (c)-(f), most II2O molecules are omitted; only those forming hydrogen bonds with the CoP 
structure are shown. Panels (e) and (f) depict the C-OH group in the trans (exo) and cis (endo) 
configuration, respectively. All snapshots are taken after 1-2 ps of short AIMD runs. Pink: Co; 
grey: C; red: O; white: H. 

For the C-OH cleavage reaction (Eq. [3]), the reaction coordinate is taken to be Rq-o^ the 
distance between the COOH carbon and the hydroxy] oxygen atom. Harmonic potentials of 
the form Eq. |8]are applied to 10 umbrella sampling windows, but with Rc-o replacing the 
4-atom coordinate R in this case. The phase space factor A{R) in the free energy expression 
analogous to Eq. [H] becomes 4:TtRq_q. 

Dealing with slow, diffusive degrees of freedom is a significant challenge in AIMD calcula- 
tions of bond-breaking free energies.—*^ During the stretching/cleavage of the C-OH bond, 



rotational phase factors of A;BTlog[47r(i?c-o)^] should naturally emerge in well-converged 
evaluations of Eq. Ei However, rotation of the nascent product 0H~ ion about the carbonyl 
carbon atom in water is too slow on AIMD time scales to accurately reproduce this rota- 
tional entropy. To give a better converged W{R), the x- and y-coordinates of the carbonyl 
C and hydroxyl O atoms are kept identical and fixed while their z coordinates are allowed 
to vary. No other atom is frozen in AIMD simulations. The rotational entropy is restored 
by multiplying AttRq^q to the Rc-o probability distribution function at each Rc_o.—^^^ 
The SI, Sec. S3, shows that this constraint has modest effect on the sampling of the Co-C-0 
angular distribution despite the bulkiness of the CoP group. 

The deprotonation free energy change in Eq. [2] is a state function which should not depend 
on the reaction coordinate chosen provided that equilibrium sampling is achieved. Using the 
one-dimensional coordinate R, however, we observe some hysteresis due to the picosecond 
time scale relaxation of the charge state of the cobalt ion as the extent of deprotonation 
varies. For Eq. [31 we are interested in the free energy barrier in addition to the free en- 
ergy change. Umbrella sampling yields a free energy barrier estimate which depends on and 
is generally underestimated by any chosen reaction coordinate because some trajectories 
with forward velocity can recross the "transition state" point and do not proceed to prod- 
uct formation.— 1^ To assess the validity of the computed barrier, we perform transmission 
coefficient (k) calculations,^'^ to be discussed in more detail below. The transition path 
sampling method^ is the rigorous approach to compute free energy barriers; although more 
costly, it will be considered in the future. Mult i- dimensional metadynamics,— new depro- 
tonation coordinates,'^ the self-consistent DFT+U method,— and new DFT functionals^^ 
may also benefit future AIMD-based electrochemical calculations, and will be mentioned in 
SeclYl 

The amount of water in the simulation cell is determined using grand canonical Monte 
Carlo simulations at constant water chemical potential, the Towhee Monte Carlo code,— 
and the SPC/E water model.— The cell size is identical to that used in AIMD simula- 
tions (13.64 A^). The temperature is set at T=300 K because the SPC/E model, unlike 
DFT/PBE, yields reasonable water structure at room temperature. One charge-neutral 
Co(III)PCOOH is placed frozen in its DFT-optimized configuration in the simulation cell. 
The CoP and COOH Lennard- Jones force field parameters are approximated with those 
of Mn(II)P— and the formate anion, respectively. The atomic partial charges are assigned 



reduced 


oxidized 


B3LYP 


PBE 


DFT-FU 


Co(I)p- 


Co(II)P 


-1.71 


-0.89 


-1.46* 


Co(I)p2- 


Co(I)p- 


-2.11 


-2.31 


NA 


Co(I)PC022" 


Co(I)PC02- 


-1.82 


-1.78* 


-1.68* 


Co(II)PCOOH- 


Co(III)PCOOH 


-1.17 


-2.14 


-1.07* 


Co(II)PCOOH2- 


Co(II)PCOOH- 


-2.17 


-2.08* 


NA 



TABLE I: Redox potential of various species using three functionals/computational methods, in 
volts. ZPE are computed using the B3LYP functional; asterisks indicate that dielectric solvation 
contributions also come from B3LYP. Convergence cannot be achieved when applying DFT-I-U to 
some doubly negatively charged systems in the gas phase. The electron affinities of these species 
are listed in the SI, Sec. S4. DFT-I-U calculations are performed using VASP; B3LYP and PBE 
calculations apply the Gaussian suite of codes. 

using Mulliken charge analysis of a gas phase B3LYP/6-311-|-G(d,p) Gaussian calculation. 
4x 10® Monte Carlo moves are attempted, 30% of them being water insertion/deletions. The 
most probable number of water molecules in the simulation cell is determined to be 71, and 
this is the water content used in AIMD simulations. A net — |e| charge and a neutralizing 
background are then imposed on the final C0PCOOH/H2O configuration from the Monte 
Carlo run, and AIMD simulations are initiated. We choose a few sampling windows along 
the reaction coordinates as seed windows. With the appropriate umbrella sampling poten- 
tials of the selected windows turned on, AIMD pre-equilibration is conducted for 2 ps at 
T=500 K. Maximally localized Wannier function analyses^ confirm that the extra electron 
resides on CoPCOOH. Then the system is further equilibrated at the target temperature 
T=425 K for another 2 ps before statistics are collected. The starting configurations for all 
other sampling windows are taken successively from snapshots of adjacent windows a few 
picosecond into their AIMD trajectories. 

The statistical uncertainty in each sampling window is estimated by splitting the trajec- 
tory into four equal parts, calculating the standard deviation for W{R2) — W{Ri), where 
Ri and R2 are the boundary values in the window, and then dividing by a/4 to yield an 
approximate error bar for the entire trajectory. The overall uncertainty convolves the stan- 
dard deviations in all windows. The statistics are generated with AIMD trajectories of at 



least 10 ps in duration; in a few windows where the uncertainties are large, 15-20 ps AIMD 
runs are conducted. 

The maximally localized Wannier function analysis^ is used to determine the charge/ spin 
state of the cobalt ion in the reactants, productions, and reaction intermediates in liquid 
water. As discussed in the SI Sec. SI, the MuUiken charge decomposition technique, pop- 
ular in the quantum chemistry community, has not been implemented in the planewave 
basis VASP code. Alternative charge decomposition methods such as simply integrating 
the charge/spin densities within some radius around the transition metal ion tend to yield 
ambiguous results.— Hence the Wannier approach appears the most useful method in our 
condensed phase, periodically replicated simulation cell setting. 



III. RESULTS 



A. Redox Potentials, Hydration Structures, and Electronic Structures 

The redox potential values are listed in Table HI To anticipate the conclusions, we find 
that the predicted B3LYP DFT+pcm absolute redox potentials are not in good agreement 
with experimental values. However, these potentials referenced to one redox couple yield 
qualitative results consistent with insights gained from AIMD simulations, and they allow 
us to assign the charge state of the key intermediate species. 

Thus, at first glance, B3LYP appears to predict a Co(I)P/Co(II)P redox couple value 
which disagrees with the average experimental value of —0.67 volt (—0.5 to —0.84, depend- 
ing on the porphyrin ring substituent and solvent^^). The redox potential also strongly 
depends on the DFT functional used, in contrast to the findings of Ref. 62] for the Co(II)P- 
NO/[Co(III)P-NO]+ couple. DFT+U and B3LYP redox potentials track each other while 
the PBE functional yields substantially different values. Systematic errors in hybrid DFT 
plus dielectric continuum redox potentials have been reported in the literature,—*^ particu- 
larly those associated with the B3LYP functional.— While part of the discrepancy in Table [T] 
may be due to DFT inaccuracies, the DFT+U redox potential, fitted to Eq. [5l is also off 
by 0.8 volt. Thus uncertainties arising from the use of the dielectric approximation used 
to calculate AGhyd as well as the lack of ring substituents in our calculations may also be 
responsible. Although PBE appears to yield a $redox for the [Co(I)P]~/Co(II)P couple close 
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FIG. 2: Mechanism of CO2 reduction with electron addition deduced from hybrid DFT plus dielec- 
tric continuum redox potential calculations. Red denotes key intermediates; green species should 
undergo fast reactions. 

to the experimental value, it performs worse for the benchmark reaction Eq. |5]than DFT+U. 
The result for that equation is deemed more reliable because the pertinent experiments were 
performed in aprotic solvents which interact weakly with the reactants.— 

To proceed, we focus on relative B3LYP redox potential values and assume that all 
electrochemical measurements occur near the B3LYP [Co (I) P]^ /Co (II) P voltage. As dis- 
cussed in the Introduction, this corresponds to the experimental condition where the onset 
of CO gas evolution is within a few tenths of a volt of the Co(II)P reduction potential. 
The B3LYP redox potentials of other species relative to [Co(I)P]^/Co(II)P (—1.71 volt) 
then determine whether an additional electron has been incorporated at that step. A sim- 



ilar approach has been used in Ref. |43|, where it has also been suggested that relative 
B3LYP redox potentials are much more reliable than absolute values. Thus the B3LYP 
functional predicts that Co(III)PCOOH is reduced to [Co(II)PCOOH]-, while [Co(I)P]- 
and [Co(II)PCOOH]- are not reduced to [Co(I)P]2- and [Co(II)PCOOH]2-, respectively, 
because their required voltages are much more negative than —1.71 volt. The experimental 
reduction potentials for [Co(I)P]~ in aprotic solvents are indeed -0.51 to -1.2 volt more neg- 
ative than that for Co(II)P.— However, unlike [Co(I)P]~ itself, [Co(I)PC02]~ is already re- 
duced to [Co(I)PC02]^~ at the [Co(I)P]"/Co(II)P voltage with a mere additional -0.11 volt. 
In other words, it is much easier to add an electron to [Co(I)PC02]~ in water than [Co(I)P]^. 



These trends can be explained with DFT+U/AIMD simulations in explicit water. The 
snapshots in Fig. [T] are generated using short, 1-2 ps AIMD trajectories, which are suffi- 
cient to yield qualitative, well-equilibrated hydration structures if not highly precise average 
hydration numbers. Despite their net negative charges, both [Co(I)P]~ and [Co(I)P]^~ are 
effectively hydrophobic plates which do not form hydrogen bonds with water molecules;— 
no water protons are observed within 2.5 A, a typical hydrogen bonding cutoff distance, of 
the porphine N and Co atoms (Figs. [T^-b). This is in contrast to the charge neutral and 
positively charged Mn(II)P and Mn(III)P, where 1 or 2 H2O molecule strongly coordinates 
to the Mn site.— CO2, a famously inert molecule, also fails to hydrogen-bond with water.— 
However, when they combine to form [Co(I)PC02]^ and [Co(I)PC02]^^, the resulting com- 
plexes form 4 to 5 hydrogen bonds with water through the partially negatively charged O 
atoms on the CO2, which now adopts a bent geometry like a carbonate^i^^i^^ or a carboxy- 
late anion2ii2^ (Figs. [U^-d). The enhanced interaction with water evidently facilitates the 
accommodation of an extra electron. This finding can be significant not just for electrochem- 
ical reduction of CO2, but also for CO2 capture in general.— The B3LYP DFT+pcm results 
reflect this COa" polarization information despite the fact that water is treated implicitly 
there. 

We also consider the hydration structures of Co(III)PCOOH and [Co(II)PCOOH]- 
(Figs. [ife-f). The C-OH trans (exo) and cis (endo) configurations are almost iso- 
energetic in gas phase Co(III)PCOOH, within 1.06 kcal/mol of each other. In contrast, 
[CoP(II)PCOOH]~ forms a strong intramolecular hydrogen bond between the COOH acid 
proton and one of the nitrogen atoms on the porphine ring when the proton is in the cis po- 
sition (Fig.[T]F). This COOH proton cannot form hydrogen bond with other water molecules. 
In the gas phase, this structure is 4.93 kcal/mol more stable than the trans configuration 
where the OH points outwards (Fig. [T]^). This feature will play a prominent role in acid-base 
reactions (Eq. [2]). Intramolecular hydrogen bonds have also been suggested to facilitate CO2 
binding and chemical reduction in the literature.-'^'^^ 

The calculated redox potentials (Table [T]) suggest an overall mechanism shown in Fig. [2l 
The CO2 adsorption and the ffist electron insertion steps are likely simultaneous and coop- 
erative. This is because the B3LYP functional predicts that the following equation, 

[Co(I)P]'- + CO2 ^ [CoP(I)PC02]'-, (9) 



exhibits a very significant free energy gain of 27 kcal/mol in water (treated as a dielectric 
continuum).— However, [Co(I)P]~ should not be readily reduced to [Co(I)P]^^ at or slightly 
below the [Co(I)P]^/Co(II)P half cell voltage, and CO2 is not strongly bound to [Co(I)P]^. 
This suggests that the CO2 may be thought of as part of the solvent, and the electron 
transfer to the cobalt complex as a solvent/C02 fluctuation-mediated process akin to the 
Marcus theory picture.— As mentioned in the introduction, the electron transfer rate may 
depend on the electrical contact between the catalyst and the gas-diffusion electrode, and is 
not the main subject of this study. The protonation and C-OH cleavage steps in Fig.[2]are at 
the heart of the catalytic function, and must be fast and spontaneous; understanding of the 
scientific principles involved can lead to improved catalysts and reaction conditions. In the 
next subsections, we consider these steps in detail using the AIMD method. Co(II)P is known 
to be weakly bound to CO.^*^ Although experimental data is not available, B3LYP/6-31+G* 
calculations indicate that the gas phase Co(I)P-CO binding energy is only 1.76 kcal/mol at 
T=0 K. Once C-OH cleavage is achieved, the rest of the reaction (Eq. H]) should proceed 
rapidly. 

The Co charge and spin states predicted for B3LYP DFT+pcm optimal structures'^ and 
in DFT+U/AIMD aqueous phase snapshots generally agree with each other. For example, 
we consider an AIMD snapshot of [Co(I)PC02]^^:H30~'" contact ion pair— in water (Fig. [3^, 
obtained using Eq. [8]with B=3 eV/A^ and Ro = —1.12 A.).^^ A maximally localized Wannier 
function analysis reveals that 6 occupied d-spin-orbitals are centered within 0.1 A of the Co 
atom; another 2 occupied spin-orbitals are 0.86 to 0.88 A away from the Co, and 1.04 and 
1.02 A away from the CO2 carbon atom, respectively. (The Co-C distance is 1.90 A in this 
snapshot.) This electronic configuration is consistent with a dative covalent bond donated 
by [Co(I)P]^^ to CO2. The singly occupied highest occupied HOMO state is delocalized on 
the porphyrin ring, as has been discussed in Paper I and is confirmed in the spin density 
plot of Fig. [3|i. We designate this species [Co(I)PC02]^^. 

A key exception to the agreement between AIMD and B3LYP DFT+pcm calculations 
is [Co(II)PCOOH]~. Figure |3]d depicts an AIMD snapshot of this species, obtained with 
harmonic constraint -Ro=— 0.4 A which yields a COOH (i.e., protonated CO2) group and 
also forces a H2O molecule to accept a hydrogen bond from the COOH proton. '^^ A Wannier 
analysis of this snapshot (Fig. |3b) reveals 7 occupied (i-spin-orbitals centered around 0.2 A 
of the Co atom. Another two occupied spin-orbitals are localized along the Co-C bond. 




FIG. 3: Spin densities as protonation of [Co(I)PC02]^ proceeds. Gold/blue: regions with net 
postive/negative spin densities. The isosurface spin density values are different for all panels to 
facilitate visualization, and are 1.32x10"^, 2.88x10"^, 2.32x10-"^, and 2.72x10"^ |e|/A3 for panels 
(a), (b), (e), and (f), respectively. The color scheme of the stick figures for water and CoP is as 
in Fig. [H (a) Well-equilibrated [Co(I)PC02]^~:H30+ contact ion pair, obtained using a harmonic 
umbrella potential (Eq. [8]) with B=3.0 eV/A, Ro = —1.12 A. The singly occupied molecular 
orbital (HOMO) state is delocalized on the porphyrin ring, (b) [Co(II)PCOOH]-, B=3.0 eV/A, 
Ro=—OA A; the harmonic potential forces a H2O molecule to accept a hydrogen bond from the 
COOH proton, (c) &: (d) reprise panels (a) and (b), respectively, omitting the spin densities 
to reveal the atomic positions more clearly. The dark blue stick figure above CoP in panel (c) 
represents the transient HsO"^. (e) & (f): Starting the trajectory from panel (a), Ro is suddenly 
switched from —1.12 A to —0.4 A. The panels are taken 1.188875 ps and 1.189125 ps into this 
trajectory and bracket the transition. 

within 0.6 and 0.8 A of the C-atom respectively. These Wannier orbital centers are closer to 
C than Co, and the electronic configuration is consistent with a [COOH]~ group attached 
to a [Co(II)P]°. The HOMO state is a Co (i-orbital; the spin density of this species (Fig. [3b) 
is far more localized than that of [Co(I)PC02]^~ (Fig[3K). In contrast, B3LYP calculations 
with dielectric approximation for water suggest a [Co(I)P]^^-COOH"'" electronic configura- 
tion, and the HOMO is a 7r-orbital on the porphyrin ring there.— This difference most likely 



reflects the explicit treatment of molecular water in the DFT+U/AIMD simulation which 
helps to stabilize a COOH" group via hydrogen bonding. In the SI, Sec. S5, explicit treat- 
ment of water molecules in the first hydration shell of the COOH group is indeed shown to 
yield a Co(II) charge state. 

The change of electronic structure coupled to protonation of [Co(I)PC02]^^ has conse- 
quences for AIMD simulations. Figures |3^-f depict the spin density transition in real time, 
occurring about 1.2 ps after the harmonic umbrella sampling potential (Eq. [8]) is suddenly 
switched from Ro = —1.12 A (consistent with that in Fig. |3K) to Ro=—OA A (Fig. [3b). As 
alluded to above, this change in Ro leads to reprotonation of [Co(I)PC02]~^ from its HsO"*" 
neighbor and the motion of a water molecule towards the C02^ group. Simultaneously, the 
electronic structure relaxes to [Co(II)PCOOH]^ due to the large driving force arising from 
the nuclear motion. However, if the driving force is not sufficient, e.g., if Ro is switched 
to an intermediate —0.7 A, the system may take much longer than 1.2 ps to spontaneously 
sample the Co(II) charge state. As discussed below, this leads to a slight hysteresis in the 
umbrella sampling calculation. 

B. Protonation of [CoPC02]^" 

Figure Hti depicts the W{R) associated with the deprotonation reaction 

[Co(I)PC02]^" + H+ ^ [Co(II)PCOOH]-. (10) 

Our 4-atom coordinate R effectively interpolates between the large negative R deprotonated 
plateau region where W{R) — ?■ 0, related to water-separated ion pairs,— and the protonated 
R > —0.4 A region where the shallow curvature of W{R) is governed by hydrogen bonding 
between the acid proton and a water molecule serving as a hydrogen bond acceptor. The dis- 
tribution of shortest distance between a water oxygen and the COOH acid proton, obtained 
in an unconstrained (i.e., i? = in the Eq. [8]t/(r)) AIMD simulation of [Co(II)PCOOH]~ in 
water, is expressed as a free energy profile {W{ro~}i)) in the inset of Fig. H^. The optimal 
ro-H is 3.2 A. Recalling the definition Eq. [7] and the fact that R2 and R^ are 0-H covalent 
bonds of lengths ~ 1 A, this optimal value implies that a true minimum in W{R) should 
not emerge until ~ 1.2 A. The optimal ro-n distance is larger than the canonical hydro- 
gen bond cutoff distance of 2.5 A, and reflects the inability of the COOH proton, engaged 



in a strong intramolecular hydrogen bond to one of the CoP nitrogen atoms, to donate a 
hydrogen bond to water molecules (Fig. [T]F). Indeed, along the AIMD trajectory, there is 
only a 2 % probability that the CO OH proton and any water oxygen atoms are within 
2.5 A of each other. Fortunately, the reaction coordinate R and the umbrella constraining 
potentials enforce hydrogen bond donation from the acid proton to water, a pre-requisite to 
deprotonation. 

To estimate pf^a via Eq. [TOl we find the most probable 0„ater-H"^ hydrogen bond dis- 
tance ro-H at each R, thus locally converting W{R) to ly(ro-H), and perform a spline fit 
to that probability distribution.—'^ The result is matched to ly(ro-H) in the small ro-H 
region obtained in the aforementioned AIMD run where the umbrella potential is absent. 
See the Fig. inset. Integrating over ro_H with a 4'7rrQ_jj volume element, which takes the 
place of the phase space factor A{R) in Eq. [6], referencing to water auto-ioniziation com- 
puted at a similar elevated temperature,— and adding a —0.57 kcal/mol zero point energy 
correction estimated from gas phase B3LYP calculations, p-ft'a=9.0±0.4 is predicted. Thus, 
[Co(II)PCOOH]^ does not behave like an ordinary carboxylate acid with pi^a ~ 4.5. The 
significant reduction of acidity indicates that protonation of [Co(I)C02]^~ is exothermic at 
the experimental pH > 7 conditions. 

In contrast, a preliminary study of the deprotonation W{R) of Co(III)PCOOII reveals 
that W{R) is almost independent of R for R < —0.4 A (not shown). Our reaction coordinate 
has been used to calculate pi^a down to 3.8 where a finite curvature persists in W{R).— 
This suggests that the pi^"a of Co(II)PCOOH is less than 3.8, much lower than that of 
[Co(II)PCOOII]~. Thus, adding an electron to the C02-ligated catalyst evidently enhances 
its ability to hold on to excess protons. While this may appear obvious in retrospect, the 
strong intramolecular hydrogen bonding in [Co(II)COOII]^ (Fig. [T]F) also likely contributes 
to its higher piCa- This preliminary p/^a estimate for Co(III)PCOOII suggests that protona- 
tion of CoPCOO^ cannot occur spontaneously at the experimental pH~7, further confirming 
the B3LYP DFT+pcm $,edox prediction that [Co(II)PCOOH]-, not [Co(III)PCOOH], is the 
key intermediate (Fig. |2]). 

Figure shows that the hydration numbers do not significantly vary with R. How- 
ever, the hydration structure may determine whether hysteresis in the cobalt charge state 
occurs as R increases. The full curve in Fig. H] collates contributions from all windows, 
such that windows 1 and 3 are initiated from [Co(I)PC02]^~ configurations generated in a 
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FIG. 4: (a) Red: Potential of mean force {W{R)) associated with [Co(II)PCOOH]~ deprotonation. 
The most negative values of R are associated with complete deprotonation while R > —0.4 A refers 
to the protonated state. Inset: the violet line depicts the W{ro^ii), computed using unbiased 
AIMD simulations; the red line depicts l^(ro-H) transformed from the biased (umbrella sampling) 
W{R) and is a cubic spline fit. (b) Hydration numbers (A^«,) of the carbonyl (green) and hydroxyl 
(blue) oxygen atoms on the COOH group. is defined as the number of water protons within 
2.5 A of these oxygen sites. Cobalt is in the Co(I) state in windows 1-3 and Co(II) state in 4- 
7, except for the dashed curves (window 3) where Co (II) prevails. The dashed lines shows that 
window 3 is affected by the cobalt charge state hysteresis (see text). 

well-equilibrated window 2 trajectory, while the window 5-7 runs originate from window 4 
([Co(II)PCOOH]^). This appears to yield a smooth W{R) curve. If the window 3 segment 
of W{R) is initiated from window 4 instead, the system remains in a Co (II) state, and the 
dashed curve, matching poorly to the window 2, materializes. This is a classic signature of 



hysteresis in umbrella sampling simulations. Fortunately, only the intermediate Ro = —0.7 A 
window 3 suffers from this problem. The maximal underestimation of [Co(II)PCOOH]~ pKa 
this can introduce is the difference between the full and dashed curves in window 3, which 
is 0.76 kcal/mol or 0.55 pH unit. This hysteresis issue is discussed in more detail in the 
appendix. 

Using a purely dielectric continuum treatment of water and the B3LYP/6-31+G* 
method, Paper I predicts deprotonation pK^, of 13.8 and 7.6 for [Co(II)PCOOH]^ and 
Co(III)PCOOH, respectively. They are several pH units higher than AIMD estimates. The 
discrepancies are partly due to the lack of explicit water molecules in the B3LYP DFT+pcm 
calculations, but they also reflect the substantial (~ 5 kcal/mol) variation in deprotonation 
free energies when using different DFT functionals (Table 4 in Paper I). 



C. C-OH bond Cleavage Reaction 

Finally, we apply umbrella sampling to study 

[Co(II)PCOOH]" ^ [Co(II)PCO] + 0H-. (11) 

We use the C-0 distance of the C-OH bond as the reaction coordinate. Figure [5^ shows 
that this reaction is exothermic. Wannier analysis reveals that the system remains in the 
Co (II) state as -Rc-o varies, and no hysteresis is observed. After accounting for standard 
state correction, the 47ri?^_o rotational contribution, a —2.5 kcal/mol ZPE correction, and 
integrating exp[—l3W{Rc-o)] in the reactant channel as in Eq. El we obtain a free energy of 
reaction AG^^^ = — 8.5±1.1 kcal/mol. The barrier height is a low AG(°)* = 5.2±0.6 kcal/mol 
confirms that the activation free energy is fairly low. We have not attempted to compute 
ZPE for AG(o)* which requires a projecton operation that removes the reaction coordinate;— 
however, in a study of a simple C-OH bond breaking reaction, ZPE has been found to be 
small, reducing AG^^^* by 0.8 kcal/mol. 
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The predicted activation barrier may depend on the -Rco coordinate chosen. As in Ref. 
we have computed the transmission coefficient.— Thus, in the umbrella sampling window 
containing the W{R) turning point, we randomly choose 10 configurations at the top of 
the barrier, half with positive velocities dRc_o/dt and half with negative ones, restart 
AIMD trajectories without umbrella sampling potentials, and determine the ratio k that 



1.5 2 2.5 3 3.5 
Rc-o (A) 



FIG. 5: (a) W(R) for the C-OH cleavage reaction, Eq. M Inset: W{R) for CO3H- ^ CO2 + 



OH from Ref. 



48l . Snapshots at points A-D are depicted in Fig. [6l (b) Hydration numbers {Ni^) 



for the carbonyl (green) and hydroxyl (blue) oxygen atoms as Rc-o varies. 

the reaction proceed without ultimate recrossing back to the reactants. k, = 1 means no 
recrossing and a perfectly chosen reaction coordinate. We find that n = 0.60, indicating 
that -Rco is a reasonable coordinate and that our reported AG* should be qualitatively 
correct. A more systematic approach is the path-sampling method, which is however more 
computationally costly. 

The C-OH cleavage activation barrier is almost a factor of 3 smaller than that previously 
found for the uncatalyzed COsH^ — )■ CO2 + 0H~ reaction (Fig. |5k inset).™ Even though 
the comparison is not perfect, in that the carbon atom is not reduced to its +2 oxidation 
state in the previous work, the cobalt porphyrin has clearly and drastically reduced the 
C-OH cleavage barrier. 

As the C-OH cleavage reaction proceeds, the carbonyl oxygen in the initially partially 
negatively charged COOH^ functional group becomes part of an uncharged carbon monoxide 



molecule (CO) weakly bound to Co(II)P, and the oxygen atom of the nascent CO exhibits a 
hydration number which steadily decreases (Fig. Eb)- In contrast, the hydroxyl oxygen 
transitions towards a hydroxide anion (OH^), and its A^^ increases to about 3.5. The A^^; 
value for the emerging 0H~ oxygen in the present heterogeneous environment is therefore 
similar to that predicted in bulk liquid water using the PBE functional. Figure [6] further 
depicts snapshots of the instantaneous hydration structure of the COOH group at different 
values of the reaction coordinate. Panel (b) represents a configuration where the COOH 
proton is intramolecularly hydrogen bonded to a porphine ring N atom (Fig. [T]F). In panel 
(a), which corresponds to a kink in Fig. [5|i, the COOH proton has instantaneously been 
donated to a N atom, forming a covalent bond with it. In the gas phase, this N-H bonded 
structure is 6.68 kcal/mol higher in energy than that of Fig. |6}3. Nevertheless, in the aqueous 
phase, this configuration is occasionally observed. 



IV. DISCUSSION 



A. Comparison with Previous CO2 Theoretical Work 

In this subsection, we make comparisons with some computational aspects of Ref. 
with Paper I.-^ 
Unlike Ref. 



48 and 



48l . we have not constrained the OH bond rotation around the C-0 axis in 
the cleaved COOH group and then corrected for the entropic contributions there. This is 
because the PBE functional we use is consistent with much faster 0H~ dynamics in water 
than the RPBE functional previously applied,— i^*^ and it is reasonable to assume that the 
OH rotation around the C-0 axis is better-sampled than in RPBE simulations within 10 ps 
AIMD trajectories. This should only affect AG^^\ not AG^^^*, because the C-OH bond is 
not completely broken at the transition state at Rc-o = 1.9 A and free OH rotation around 
the C-0 axis does not occur there. The important qualitative conclusion of this work is 
that AC^^^ is exothermic and that C-OH bond breaking is thermodynamically downhill; the 
precise free energy change associated with Eq. [11] is less important. 

We have not attempted to correct the AIMD/DFT-FlJ W{R) with single point MP2 



calculations as was done in Ref. 
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The CoP systems examined in this work are too large 
for large-basis MP2. Furthermore, in Ref. ^s], the DFT functional used was RPBEj^SS which 




(d) 



FIG. 6: (a)-(d): snapshots in the four umbrella sampling windows indicated in Fig. [5j AIMD 
simulations are conducted in an explicit liquid water environment; only a few water molecules are 
shown in these snapshots. In panel (a), the proton on the COOH group has migrated to one of the 
nitrogen atoms on the porphyrin ring. The color scheme is as in Fig. [TJ 

is arguably less accurate for heterogeneous C-0 bond breaking than the PBE functional 
used herein. 

Unlike Paper I, we have focused on Eq. [3] and not the proton-assisted variation 



Thermodynamically the two are equivalent.— In terms of kinetics, which Paper I does not 
deal with, they will have different activation barriers. To estimate AG'-'^^* for Eq. [121 we note 
that all AIMD W{R) reported for protonation/deprotonation reactions in the literature have 
been monotonic, i.e., the free energy changes and activation barriers are the same. Assuming 
the second half of Eq. [12] is also fast to the point of being barrierless, which is a lower limit 
on the overall Eq. [12] barrier, AG'^^^* for Eq. [4] would be entirely due to the AG^^^ of the first 
half of this equation and related to the pKa of Co(II)PCO(OH2). At the pll=7 experimental 
conditions, the pKg, of the C-OH2 group in Co(II)PCO-OH2 will have to be above 3.2 in 
order to have a lower barrier than the AG* = 5.2 kcal/mol we find for Eq. [3] This appears 
unlikely; C-OH2 groups tend to be very acidic and lose protons readily. In any case, if this 



[Co(II)PCOOH]" + H+ ^ Co(II)PCO(OH2); 
Co(II)PCO(OH2) ^ [Co(II)PCO] + H2O. 



(12) 



alternate, proton assisted route of C-0 cleavage were faster, the barrier estimated in our 
PMF calculation (Fig. E]) would be an upper bound to the reaction activation free energy; 
the C-0 bond breaking step would still exhibit fast dynamics, and the qualitative conclusion 
of this paper would be unchanged. We plan to revisit Eq. [12] in the future. 

B. Alternative Computational Methods 

A two-dimensional PMF calculation should remove the hysteresis behavior in the pro- 
tonation reaction (Eq. [Ill Fig- HI)- A convenient second variable may be the hydration 
number of the O atom in the OH group. This is because the black dashed curve in Fig. 
exhibits A^=1.49 for that oxygen, considerably higher than the A^=0.75 for the red, well- 
behaved W{R) segment, suggesting that the electronic structure is correlated with the aver- 
age hydration number. 2-D PMF simulations would be substantially accelerated using the 
metadynamics method^^"^^ over traditional umbrella sampling. 

Our pf^a calculations suggest that deprotonation of weak acid groups which exhibit in- 
tramolecular hydrogen bonding and do not donate hydrogen bond to water may require 
the use of a large number of sampling windows. This difficulty may be circumvented by 
reversibly annihilating the proton using an artificial reaction pathway^^ in a what might be 
called a "molecular grand canonical Monte Carlo" approach.— This method does introduce 
the disadvantage of changing the net charge in the finite-sized simulation cell, and may 
require a number of new conformational constraints. 

Finally, the self-consistent DFT+U approach^^ has been tested for Eq. [5] This method has 
the potential to establish a U value without resorting to parameterization with experimental 
results. Our preliminary studies suggest that this approach yields a Co(III)P-CO binding 
energy that is too small, but further development of this promising approach is under way— 

The above theoretical considerations, touching on many newly developed techniques, 
emphasize the complexity and challenges associated with modeling electrochemical reactions 
in explicit-water aqueous phase simulations. 



C. Accuracy of DFT functionals and redox potentials 



Accurate DFT functionals and dielectric continuum approximation of the aqueous solvent 
are critical for DFT+pcm determination of redox potentials,— which in turn govern the 
viable reaction intermediates and the overall reaction mechanism. We have so far considered 
PBE, B3LYP, and DFT+U electronic structure methods with the U parameter in DFT+U 
fitted to an experimental binding constant (Table [T]). B3LYP and DFT+U redox potentials 
track each other and should predict the same reaction mechanism at voltages slightly more 
negative than the Co(I)P/Co(II)P couple. In contrast, the PBE functional predicts that 
[Co(II)COOH]^ is extremely unstable with respect to the far more acidic [Co(III)COOH] 
complex near the PBE Co(I)P/Co(II)P redox potential. The B3LYP $redox predictions 
are more consistent with AIMD hydration structure considerations and the fact that the 
C02-reduction reaction readily occurs near neutral pH in a carbonate buffer (Sec. IIIB . 

More accurate DFT functionals may be used in the future to further examine the mech- 
anism proposed herein. Candidates include the M06 class of functionals designed to yield 
better thermochemistry accuracy for transition metal complexes,— the B4(XQ3)LYP func- 
tional which has been found to yield improved redox potentials for a suite of test cases,— 
and Gutzwiller wavefunction based methods.— The quality of the PCM dielectric continuum 
model2^ should also be further examined.— 

V. CONCLUSIONS 

In this work, we have applied first principles calculations to examine the mechanism of 
the multi-step, two-electron electrochemical reduction of CO2 to CO in water using cobalt 
porphyrin (CoP) as catalyst. First we have extracted redox potentials from DFT plus dielec- 
tric continuum solvation calculations using the B3LYP functional and a dielectric continuum 
treatment of water.— Even though the absolute value for the [Co(I)P]^/Co(II)P couple is not 
in good agreement with experiments, the relative values of various redox potentials allow us 
to determine where the electron transfers occur among the four intermediate steps. Due to 
the enhanced interaction of CO2 with water when bound to cobalt porphine, [Co(I)PC02]^~ 
and [Co(II)PCOOII]^ are the key intermediates. This finding may be useful not just for 
electrochemical reduction of CO2, but for CO2 capture from flue gas as well. 



AIMD umbrella sampling calculations show that the pK^ associated with 
[Co(II)PCOOH]~ deprotonation is about 9.0. This indicates that the protonation of 
[Co(l)PC02]^~ is downhill at the bicarbonate buffer experimental conditions (pH ~7). The 
subsequent cleavage of the C-OH bond is also exothermic, and the activation free energy 
involved is estimated to be only 5.2 kcal/mol. If we assume a vibrational pre-factor of 
k=0.1 ps~^, C-OH cleavage should occur in nanosecond timescale at T=300 K. Hence two 
key steps in the multistep reaction should proceed readily, and it is likely that the electron 
transfer between the gas diffusion electrode and the polymerized porphyrin catalyst is the 
rate limiting step of the CO2 to CO reduction reaction in water.-"— 
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Appendix: Hysteresis in the Deprotonation AIMD Simulation 

This appendix discusses in more detail the hysteresis behavior during pKa PMF calcula- 
tions (window 3 in Fig. HJ. 

As mentioned in Sec. Ull all AIMD trajectories ultimately originate from a Monte 
Carlo simulation-equilibrated classical force field configuration where a CoPCOOH is fixed 
in a gas phase optimized, intramolecularly hydrogen-bonded geometry (Fig. [IT). The 
[Co(II)PCOOH]^ is immersed in water and equilibrated using AIMD at several stretched 
values of the reaction coordinates associated with Eq. and Eq. [H] using harmonic poten- 
tials U{R) (Eq.E]). 

During the equilibration run for the deprotonated window 2 of Fig. HI we impose Ro = 
-1.12 A on the initial [Co(II)PCOOH]- system, which turns into a [Co(I)PC02]2":H30+ 
contact ion pair within 1 ps. A maximally localized Wannier function analysis of an AIMD 
snapshot confirms that, in this ion pair, Co has spontaneously switched to the +1 charge state 
which should be favored for large negative R. Therefore the AIMD trajectory in window 2 
(and that in window 1, spawned from window 2) yields an unambiguous Co(I) charge state. 
Likewise, windows 4-7 are spawned successively from a window 4 trajectory equilibrated 
using Ro = —0.4 A, and they reflect a Co(II) charge state which is favorable in this R range. 
An additional test further confirms that the Co (II) charge state spontaneously occurs at less 
negative R. We start with a [Co(I)C02]^~, Ro = —1.12 A configuration in a equilibrated 
window 2 trajectory and abruptly switch to Ro = —0.4 A (i.e., effectively jumping from 
window 2 to window 4 in Fig. H]). [Co(II)PCOOH]~ is recovered within 1.2 ps, consistent 
with the Co(II) charge state observed in the window 4 trajectory which originally started 
out as [Co(II)PCOOH]~. The change in electronic structure in real time, represented by 
the changes in the spatial distribution of the spin density, is depicted in Fig. [3l Thus, after 
an approximately 1 ps equilibration run, the final Co charge state becomes independent of 
initial conditions. 

Only window 3, located in the Co(I)/Co(II) transition region, exhibits a strong depen- 
dence on initial conditions. We have initially started from a snapshot in window 4 and 
then switched Ro = —0.4 A to -Ro = —0.7 A. The latter value of Ro does not apparently 
contain sufficient driving force to rapidly alter the Co charge state, and the system re- 
mains [Co(II)PCOOH]~ throughout the 10 ps sampling trajectory. This yields the dashed 



W{R) segment in Fig. which exhibits a slope at the window edge that matches poorly 
to the window 2, [Co(I)PC02]^~:H30+ contribution. To make progress, we restart the win- 
dow 3 simulation from a snapshot of window 2, abruptly switch the umbrella potential from 
Ro = —1.12 A there to Ro=-0.7 A, equilibrate for 1 ps, and collect statistics for 10 ps. 
The system remains in the Co(I) charge state throughout the trajectory. The corresponding 
W{R) segment in window 3 (full curve in Fig. H^) matches reasonably well with those in 
windows 2 and 4, and is taken to be the final result. 

An ergodic AIMD simulation should in principle spontaneously sample both Co (I) and 
Co(II) charge states. Thus, the correct W{R) in this intermediate R region should be a 
weighted average of the full and dashed curves. Apparently the statistical weight for Co(I) 
is much larger, so that only including Co(I) W{R) information already yields a smooth 
W{R) curve in Fig. HI 

A secondary, less signficant type of hysteresis associated the position of the COOH proton 
also becomes apparent in Fig. [71 When the COOH group is intact, and no H2O accepts 
a hydrogen bond from the COOH proton, the hydroxyl group intramolecularly hydrogen- 
bonds to a N atom on the porphine ring in a cis configuration (Fig. [If). Enforcing U{R) with 
R < 0.5 A imposes COOH-H2O hydrogen bonding that breaks this intramolecular coupling 
(Fig. [7^)- When R is further reduced to i? ~ — 1 A, the [COOH]^ proton is detached 
from the hydroxyl oxygen, landing on the hydrogen-bonding accepting H2O molecule. This 
newly formed HsO^ spontaneously migrates away from the hydrophobic porphine ring and 
coordinates to what is now a COg" group in the axial (trans) position. Now the excess 
proton sticks out of the porphine plane. When we increase Ro to reprotonate the COOH 
group, the system remains in this isomeric form (Fig. Wp)- Isomerization between the two 
may entail a free energy barrier of several kcal/mol even for a stretched 0-H, and does not 
occur on AIMD timescales in window 3. The lack of isomerization should have no significant 
effect on the W{R) of window 3, while the gas phase Fig. [If [Co(II)PCOOH]^ intramolecular 
hydrogen bonded configuration is stabilized over the Fig. [1^ isomer by about 4.93 kcal/mol, 
that hydrogen bond is already broken when the carboxylate proton donates a hydrogen bond 
to a water molecular (Fig. [7^). This secondary hysteresis may be avoided altogether using 
an artificial deprotonation coordinate,— although it is not obvious the electronic hysteresis 
will be avoided as well. 



FIG. 7: Panels (a) & (b) correspond to the solid and dashed line in window 3, respectively. The 
color scheme is as in Fig. [TJ 
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